// Copy of yyt.h - 28.03.2011 - using Mathematica

/*
 *********************************************************
 ******************** Structures *************************
 *********************************************************
*/

struct parameters {
	double sdelta;
	double J;
	double t;
	double U;
	unsigned int L;
};

struct variables {
	double BDelta;
	double chi;
	double m;
	double d;
	double mu;
};

/*
*********************************************************
***************** Small useful functions ****************
*********************************************************
*/

void clear(parameters &p) {
	p.sdelta = p.J = p.t = p.J = p.L = 0;
}

void clear(variables &v) {
	v.BDelta = v.chi = v.m = v.d = v.mu = 0;
}

void addTo(variables &v, variables &q, double d=1.) {
	q.BDelta = v.BDelta*d;
	q.chi = v.chi*d;
	q.m = v.m*d;
	q.d = v.d*d;
	q.mu = v.mu*d;
}

inline double abs(double a) {
	if(a<0) a=-a;
	return a;
}

/*
 *********************************************************
 ***************** CLASS Functions ***********************
 *********************************************************
*/

class Func {
private:
	parameters p; // parameters
	variables v; // variables
	double n, r, w; // contractions

	void refresh() {
		n = 1.-p.sdelta;
		r = n/2. + v.m;
		w = n/2. - v.m;
	}

public:
	Func(const parameters &p, const variables &v) : p(p), v(v) { refresh(); }

	void setPar(const parameters &p) { this->p = p; refresh(); }
	void setVar(const variables &v) { this->v = v; refresh(); }

	double gt();
	double gs();
	double dgtdm();
	double dgtdd();
	double dgsdm();
	double dgsdd();
	double gamma(double kx, double ky);
	double eta(double kx, double ky);
	double epsilon(double kx, double ky);
	double BDelta_d();
	double BDelta_af();
};

double Func::gt() {
	return ((sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.)) +
		       sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.)))*
		     (sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.)) +
		       sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.)))*(-2*v.d + n))/(-2*(-v.m + n/2.)*(v.m + n/2.) + n);
}

double Func::gs() {
	return pow(-2*v.d + n,2)/pow(-2*(-v.m + n/2.)*(v.m + n/2.) + n,2);
}

double Func::dgtdm() {
	return -(((sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.)) +
	          sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.)))*(-2*(-v.m + n/2.) + 2*(v.m + n/2.))*
	        (sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.)) +
	          sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.)))*(-2*v.d + n))/
	      pow(-2*(-v.m + n/2.)*(v.m + n/2.) + n,2)) +
	   ((sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.)) +
	        sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.)))*
	      ((-(((1 + v.d - n)*(1 - v.m - n/2.))/pow(1 + v.m - n/2.,2)) -
	           (1 + v.d - n)/(1 + v.m - n/2.))/
	         (2.*sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.))) +
	        (v.d/(-v.m + n/2.) + (v.d*(v.m + n/2.))/pow(-v.m + n/2.,2))/
	         (2.*sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.))))*(-2*v.d + n))/
	    (-2*(-v.m + n/2.)*(v.m + n/2.) + n) +
	   ((((1 + v.d - n)/(1 - v.m - n/2.) +
	           ((1 + v.d - n)*(1 + v.m - n/2.))/pow(1 - v.m - n/2.,2))/
	         (2.*sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.))) +
	        (-((v.d*(-v.m + n/2.))/pow(v.m + n/2.,2)) - v.d/(v.m + n/2.))/
	         (2.*sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.))))*
	      (sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.)) +
	        sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.)))*(-2*v.d + n))/(-2*(-v.m + n/2.)*(v.m + n/2.) + n);
}

double Func::dgtdd() {
	return (-2*(sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.)) +
	        sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.)))*
	      (sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.)) +
	        sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.))))/(-2*(-v.m + n/2.)*(v.m + n/2.) + n) +
	   ((sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.)) +
	        sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.)))*
	      ((1 - v.m - n/2.)/
	         (2.*sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.))*(1 + v.m - n/2.)) +
	        (v.m + n/2.)/(2.*(-v.m + n/2.)*sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.))))*(-2*v.d + n))/
	    (-2*(-v.m + n/2.)*(v.m + n/2.) + n) +
	   (((1 + v.m - n/2.)/
	         (2.*(1 - v.m - n/2.)*sqrt(((1 + v.d - n)*(1 + v.m - n/2.))/(1 - v.m - n/2.))) +
	        (-v.m + n/2.)/(2.*sqrt((v.d*(-v.m + n/2.))/(v.m + n/2.))*(v.m + n/2.)))*
	      (sqrt(((1 + v.d - n)*(1 - v.m - n/2.))/(1 + v.m - n/2.)) +
	        sqrt((v.d*(v.m + n/2.))/(-v.m + n/2.)))*(-2*v.d + n))/(-2*(-v.m + n/2.)*(v.m + n/2.) + n);
}

double Func::dgsdm() {
	return (-2*(-2*(-v.m + n/2.) + 2*(v.m + n/2.))*pow(-2*v.d + n,2))/
			   pow(-2*(-v.m + n/2.)*(v.m + n/2.) + n,3);
}

double Func::dgsdd() {
	return (-4*(-2*v.d + n))/pow(-2*(-v.m + n/2.)*(v.m + n/2.) + n,2);
}


double Func::gamma(double kx, double ky) {
   return 2.*(cos(kx) + cos(ky));
}

double Func::eta(double kx, double ky) {
   return 2.*(cos(kx) - cos(ky));
}

double Func::epsilon(double kx, double ky) {
   return -(gt()*p.t + 3./8.*gs()*p.J*v.chi)*gamma(kx,ky);
}

double Func::BDelta_d() {
   return 3./8.*gs()*p.J*v.BDelta;
}

double Func::BDelta_af() {
   return 2.*gs()*p.J*v.m;
}

/*
 *********************************************************
 ******************* CLASS Equations *********************
 *********************************************************
*/

class Equa {
private:
	parameters p; // parameters
	variables v; // variables
	Func *f;

	// result of f (independent on k_x and k_y)
	double gt;
	double gs;
	double dgtdm;
	double dgtdd;
	double dgsdm;
	double dgsdd;
	double BDelta_d;
	double BDelta_af;

	void onePoint(variables &y, int i, int j);

	void refreshRF() {
		gt = f->gt();
		gs = f->gs();
		dgtdm = f->dgtdm();
		dgtdd = f->dgtdd();
		dgsdm = f->dgsdm();
		dgsdd = f->dgsdd();
		BDelta_d = f->BDelta_d();
		BDelta_af = f->BDelta_af();
	}

public:
	Equa(parameters &p, variables &v) : p(p), v(v) { f = new Func(p, v); refreshRF(); }
	~Equa() { delete f; }

//	void setPar(parameters &p) {
//		this->p = p;
//		f->setPar(p);
//		refreshRF();
//	}

	void run(variables &w);
};

void Equa::onePoint(variables &w, int i, int j) {
	double x = (double)i*2.*M_PI/p.L;
	double y = (double)j*2.*M_PI/p.L;

	double gamma = f->gamma(x, y);
	double eta = f->eta(x, y);
	double epsilon = f->epsilon(x, y);
	double dzeta = sqrt(epsilon*epsilon + BDelta_af*BDelta_af);
	double a1 = dzeta-v.mu;
	double a2 = -dzeta-v.mu;
	double b = BDelta_d*eta;
	double E1 = sqrt(a1*a1 + b*b);
	double E2 = sqrt(a2*a2 + b*b);

		double dE1dgt = -((pow(gamma,2)*p.t*((-3*v.chi*p.J*gs)/8. - gt*p.t)*
			       (-v.mu + sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) +
			           pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2))))/
			     (sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) + pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2))*
			       sqrt((9*pow(v.BDelta,2)*pow(eta,2)*pow(p.J,2)*pow(gs,2))/64. +
			         pow(-v.mu + sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) +
			             pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2)),2))));

		double dE1dgs = ((9*pow(v.BDelta,2)*pow(eta,2)*pow(p.J,2)*gs)/32. +
			     ((8*pow(p.J,2)*pow(v.m,2)*gs - (3*v.chi*pow(gamma,2)*p.J*((-3*v.chi*p.J*gs)/8. - gt*p.t))/4.)*
			        (-v.mu + sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) +
			            pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2))))/
			      sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) + pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2)))/
			   (2.*sqrt((9*pow(v.BDelta,2)*pow(eta,2)*pow(p.J,2)*pow(gs,2))/64. +
			       pow(-v.mu + sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) +
			           pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2)),2)));

		double dE2dgt = (pow(gamma,2)*p.t*((-3*v.chi*p.J*gs)/8. - gt*p.t)*
			     (-v.mu - sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) +
			         pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2))))/
			   (sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) + pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2))*
			     sqrt((9*pow(v.BDelta,2)*pow(eta,2)*pow(p.J,2)*pow(gs,2))/64. +
			       pow(-v.mu - sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) +
			           pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2)),2)));

		double dE2dgs = ((9*pow(v.BDelta,2)*pow(eta,2)*pow(p.J,2)*gs)/32. -
			     ((8*pow(p.J,2)*pow(v.m,2)*gs - (3*v.chi*pow(gamma,2)*p.J*((-3*v.chi*p.J*gs)/8. - gt*p.t))/4.)*
			        (-v.mu - sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) +
			            pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2))))/
			      sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) + pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2)))/
			   (2.*sqrt((9*pow(v.BDelta,2)*pow(eta,2)*pow(p.J,2)*pow(gs,2))/64. +
			       pow(-v.mu - sqrt(4*pow(p.J,2)*pow(v.m,2)*pow(gs,2) +
			           pow(gamma,2)*pow((-3*v.chi*p.J*gs)/8. - gt*p.t,2)),2)));

	double part = (dzeta-v.mu)/E1 + (dzeta+v.mu)/E2;
	double dEvardt = -dE1dgt - dE2dgt;
	double dEvards = -dE1dgs - dE2dgs;

	w.BDelta += eta*eta*BDelta_d*(1./E1 + 1./E2);
	w.chi += gamma*epsilon/dzeta*((-dzeta+v.mu)/E1 + (-dzeta-v.mu)/E2);
	w.m += BDelta_af/dzeta*part - (dEvardt*dgtdm + dEvards*dgsdm)/(2.*gs*p.J);
	w.d += dEvardt*dgtdd + dEvards*dgsdd;
	w.mu += (dzeta-v.mu)/E1 + (-dzeta-v.mu)/E2;
}

void Equa::run(variables &w)
{
	clear(w);
	variables temp; clear(temp); // temporary vector using during summation

	unsigned int max = p.L/2; // attention: L should be even!

	// SMART WAY OF SUMMATION - I use symmetry k_x,y <-> -k_x,y
	// * edges (2 times)
	for(unsigned int i=1; i<max; i++) {
		onePoint(temp, i, 0);
		onePoint(temp, 0, i);
		onePoint(temp, i, max-i);
	}
	addTo(temp, w, 2);
	clear(temp);
	// * interior area (4 times)
	for(unsigned int i=1; i<max; i++) {
		for(unsigned int j=1; j<i; j++) {
            onePoint(temp, max-i, j);
		}
	}
	addTo(temp, w, 4);
	clear(temp);
	// * missing corners
	onePoint(w, max, 0);
	onePoint(w, 0, max);
	onePoint(w, 0, 0);

	double chi2 = v.chi*v.chi;
	double m2 = v.m*v.m;
	double BDelta2 = v.BDelta*v.BDelta;

	double L2 = (double)(p.L*p.L);
	w.BDelta /= 4.*L2;
	w.chi /= 4.*L2;
	w.m /= 2.*L2;
	w.d /= L2;
	w.mu /= L2;

	w.BDelta += -v.BDelta;
	w.chi += -v.chi;
	w.m += -v.m - dgsdm*(3./4.*(BDelta2+chi2) + 2.*m2)/(4.*gs);
	w.d += p.U + dgsdd*(3./4.*(BDelta2+chi2) + 2.*m2)*p.J;
	w.mu += -p.sdelta;

	if(w.BDelta != w.BDelta) { w.BDelta = 1000; }
	if(w.chi != w.chi) { w.chi = 1000; }
	if(w.m != w.m) { w.m = 1000; }
	if(w.d != w.d) { w.d = 1000; }
	if(w.mu != w.mu) { w.mu = 1000; }
}
